#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <sys/time.h>
#include <cuda.h>


/*
 *  Function: SmoothField
 *
 *  Smooth a field of dimension dim, dim using a kernel of width width
 *  pFieldIn should point to an input array of size dim*dim
 *  pFieldOut should point to an output array of size dim*dim.
 *     One should initialize(memset) pFieldOut to 0 so that the border
 *     do not contain non-zero values
 */
__global__ void SmoothField ( unsigned dim, unsigned width, float* pFieldIn, float* pFieldOut )
{

	unsigned int row = blockIdx.y * blockDim.y + threadIdx.y;
	unsigned int col = blockIdx.x * blockDim.x + threadIdx.x;
	//Go over kernel and get total value and average
	unsigned int window_width = (2*width+1);
	if(i < dim && j < dim)
	{
		/**
		float value = 0;
		int num_kernel_cells = window_width * window_width;	
  		for ( unsigned int i=0; i<window_width; i++ )
  		{
			int kernel_row = (row + i) * dim;
    			for ( unsigned int j=0; j<window_width; j++ )
    			{
				int kernel_col = (col + j);
      				value += pFieldIn[kernel_row + kernel_col];
    			}
  		}
		pFieldOut[*dim + col] = blockIdx.y*blockDim.y;
		**/
		pFieldOut[i*dim + j]
	}
}

/*
 * Function: InitField
 *
 *  Initiliaze the input field to well known values for debugging
 *  and visualization purposes.
 */
__host__ void InitField ( unsigned dim, float* pField )
{
  for ( unsigned i=0; i<dim; i++ )
  {
    for ( unsigned j=0; j<dim; j++ )
    {
      pField[i*dim+j] = j + i;
    }
  }
}

/*
 * Function: main
 *
 *  See inline comments describing steps.
 */
int main ()
{
  // Set the parameters of the kernel smoother.
  //  the valid region is from
  //   corner (kernelwidth, kernelwidth) to
  //   corner (dimension-kernelwidth-1, dimension-kernelwidth-1)
  unsigned dimension = 2064;
  unsigned kernelwidth = 8;
  unsigned smoothed_dim = (dimension-(2*kernelwidth));
  unsigned num_cells = smoothed_dim * smoothed_dim;

  float *a_new, *a_prev, *a_prev_host, *a_new_host; //Pointer to new and previous matrix

    size_t size_prev = dimension * dimension * sizeof(float);
    cudaMalloc((void **) &a_prev, size_prev);   // Allocate array on device
    a_prev_host = (float *)malloc(size_prev);

  InitField(dimension, a_prev_host);
  cudaMemcpy(a_prev, a_prev_host, size_prev, cudaMemcpyHostToDevice );
/**  
for ( unsigned i=0; i< dimension; i++ )
  {
    for ( unsigned j=0; j< dimension; j++ )
    {
      printf ("%d %d %f\n", i, j, a_prev_host[i*dimension + j]);
    }
 }
**/  
	size_t size_new = num_cells * sizeof(float);
    	cudaMalloc((void **) &a_new, size_new);   // Allocate array on device
	a_new_host = (float *)malloc(size_new);
  
  	cudaMemset ( a_new, 0, num_cells *sizeof(float));
  	int block_size = 16;
	dim3 dimBlock(block_size,block_size);
	dim3 dimGrid((smoothed_dim+dimBlock.x-1)/dimBlock.x, (smoothed_dim+dimBlock.y-1)/dimBlock.y);
  // Collect timing information
  struct timeval ta, tb;
  gettimeofday ( &ta, NULL );

  // Invoke the kernel smoother
	printf("hello\n");
  SmoothField <<< dimGrid, dimBlock>>> ( smoothed_dim, kernelwidth, a_prev, a_new );
cudaThreadSynchronize();
	printf("goodbye\n");
  // Report timing information
  gettimeofday ( &tb, NULL );

  if ( ta.tv_usec < tb.tv_usec )
  {
    printf ("Elapsed total time (s/m): %li:%d\n", tb.tv_sec - ta.tv_sec, tb.tv_usec - ta.tv_usec );
  } else {
    printf ("Elapsed total time (s/m): %li:%d\n", tb.tv_sec - ta.tv_sec - 1, 1000000 - tb.tv_usec + ta.tv_usec );
  }

  cudaMemcpy (a_new_host, a_new, size_new,cudaMemcpyDeviceToHost);
  //See what happened.

  for ( unsigned i=0; i < smoothed_dim; i++ )
  {
    for ( unsigned j=0; j < smoothed_dim; j++ )
    {
      printf ("%i %i %f \n",i,j,a_new_host[(i*smoothed_dim) + i]);
    }
  }
	printf("dimBlock.x=%d dimGrid.x=%d\n",dimBlock.x,dimGrid.x);
	printf("dimBlock.y=%d dimGrid.y=%d\n",dimBlock.y,dimGrid.y);
  cudaFree(a_new);
  cudaFree(a_prev);
  free(a_prev_host);
  free(a_new_host);
  return 0;	
}
